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Motivated by recent work on instabilities in expanding domains in reaction-diffusion settings, 
we propose an analog of such mechanisms in energy-conserving wave equations. In particular, 
we consider a nonlinear Schrodinger equation in a finite domain and show how the expansion or 
contraction of the domain, under appropriate conditions, can destabilize its originally stable solutions 
OO ■ through the modulational instability mechanism. Using both real and Fourier space diagnostics, 

' we monitor and control the crossing of the instability threshold and, hence, the activation of the 

, instability. We also consider how the manifestation of this mechanism is modified in a spatially 

CN| . inhomogeneous setting, namely in the presence of an external parabolic potential, which is relevant 

^ ' to trapped Bose-Einstein condensates. 
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I. INTRODUCTION 

There has been a considerable recent interest in the study of pattern formation in expanding domains, and pattern 
formation supported by them, in the context of reaction-diffusion equations (see, e.g., [l|,[2|,S01 and references therein 
for a substantial body of earlier work). Such mechanisms are viewed in mathematical biology as a possibly important 
element in the formation of skin patterns on species of animals or fishes [5]. One of the key ideas in that regard is 
C ■ that in the case of finite-size domains, Turing- type mechanisms [6] are only activated for appropriate wavenumbers, 
which, in turn, signify appropriate sizes of the domain. In particular, as is pointed out in [1], the domain size should 
be such so as to admit at least a half-cycle of the intrinsic pattern (to be formed by the instability) wavelength. Once 
activated, these mechanisms typically lead to the destabihzation of a uniform concentration profile, due to the presence 
(S| ■ of diffusion, and the formation of localized waveforms through frequency-doubling and related mechanisms [1, 7]. 
^ ■ Our aim in the present work is to illustrate that such mechanisms are not only relevant in reaction-diffusion 
^ : equations and pattern- forming dissipative systems, but also to Hamiltonian nonlinear dispersive wave equations, such 
as the prototypical model of the nonlinear Schrodinger (NLS) equation [8]. The latter class of models also features 
instabilities which emerge only for appropriate bands of wavenumbers; a prominent example of such instabilities is 
the modulational instability (MI) of the plane wave solutions of the focusing NLS equation (see e.g. ^] for a recent 
I review and for a detailed discussion). This instability only emerges (for a fixed amplitude of the plane wave) for 
sufficiently small wavenumbers, i.e., for sufficiently large perturbation wavelengths, or (for a fixed wavenumber) for 
I sufficiently large amplitudes of the plane wave. Hence, one can envision a setting where the initial domain is small 
enough that the smallest wavenumber permissible by the domain (say, 27r/L for periodic boundary conditions and 
domain size L) is larger than the critical wavenumber for the onset of the instability. Equivalently, the domain size 
may be smaller, at time t = 0, than the minimum wavelength necessary for the instability. Then, with an appropriate 
■ temporal dependence, the domain size may increase and cross the relevant instability threshold; this, in turn, leads 
5^ I to the emergence of patterns in a way reminiscent of that occur ing in their Turing analogs in dissipative systems. In 
fact, these patterns in the Hamiltonian case will consist of robust solitary waves, as we will illustrate below, since MI 
is one of the key mechanisms leading to the formation of such structures 0, ^3] • 

This mechanism could be of relevance to a wide array of Hamiltonian systems, since the NLS model is relevant to 
a diverse set of applications including the propagation of hght in optical fibers [l^ , the dynamics of Bose-Einstein 
condensates (BECs) at extremely low temperatures [ijj, plasma physics and fluid dynamics [l2|], and so on. In fact, 
MI in itself is still one of the most active topics of investigation in some of these contexts. In particular, a variant 
of MI was used in the first experiment that demonstrated the generation of matter- wave soliton trains in BECs [l^ . 
Furthermore, MI is recently becoming increasingly accessible in a variety of settings experimentally (and theoretically); 
these include spatially periodic ones such as BECs confined in optical lattices (where MI introduces a classical form 
of a superfluid-insulator transition [13]) and optical waveguide arrays [15], as well as ones that are periodic in the 
evolution variable, such as the nonlinear layered optical media recently studied in Refs. [16]. 

The presentation of our results will be structured as follows. We will first give a brief theoretical background on 
the nature of the instability and explain how it is affected by the existence of a finite domain size. We will then 
proceed to examine four different numerical experiments, attempting to elucidate the emergence of the instability for 
time-dependent domain sizes. Our first and most straightforward example will be the one discussed above: a stable 
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uniform solution (for sufficiently small domain size) will be subjected to a time-dependent domain expansion, resulting 
in what we will call a "soliton sprinkler" through the emergence of ML Our second experiment will be somewhat 
counter-intuitive in that the instability will emerge even though the domain size is decreasing. This will be achieved 
by preserving the norm of the solution (its power in optics or its number of atoms in BEG), and taking advantage 
of the instability dependence on the beam intensity. The last two examples will introduce an expanding or contracting 
domain situation in a somewhat modified but also directly relevant and more straightforwardly realizable experimental 
setting, inspired by the physics of BECs. In particular, we will use a parabolic potential (routinely employed to trap 
the BEG) with a time-dependent frequency, which, in turn, may weaken or strengthen the confinement of the atoms, 
introducing a corresponding effective expansion or contraction of the domain. We will seek signatures of MI in this 
more complex, yet more realistic setting. Finally, we will summarize our findings and present our conclusions. 



II. RESULTS 



A. Theoretical Background 

We consider the following dimensionless NLS equation for the field u{x^t), 

idtu = -dlu + g\u\^u + l/ext(^)^, (1) 

where g is the nonlinearity coefficient and lext(^) is an external potential. Note that in the context of BEGs this equa- 
tion is usually referred to as the Gross-Pitaevskii (GP) equation [lH); in that case, u{x^t) is the BEG wavefunction, 
g is proportional to the s-wave scattering length, and Ve^t{x) is a trapping potential where the BEG is confined. For 
simplicity, the nonlinearity coefficient is scaled so that ^ = ±1, corresponding, respectively, to repulsive or attractive 
interatomic interactions (or, generally, to defocusing or focusing nonlinearity). Equation ^ conserves the energy, as 
well as the squared norm, N = \u\'^dx (in the context of BEGs, N represents the normalized number of atoms 
of the condensate). 

For our theoretical analysis, we first consider the homogeneous (untrapped) case, i.e., Eq. ^ with Ve^t = 0, 
which admits plane-wave solutions of the form u{x^t) = uoex.p[i{kx — ut)], where uq^ /c, and u are, respectively, 
the amplitude, wavenumber and frequency of the plane wave solution. These parameters are connected through the 
dispersion relation: cj = k'^-\-guQ. The stability of the aforementioned solution is analyzed by introducing the following 
linear stability ansatz into Eq. ([T]): 

u{x^ t) = {uq + eb) ex.-p[i{{kx — uot) + ea)], (2) 

where h and a are perturbations of the amplitude and the phase respectively. By solving the resulting linear equations 
[to 0(e)] through the Fourier mode expansion a = aoexp[z((5a; — 77^)] and h = 6oexp[z((5x — 77^)], where Q and rj are 
the wavenumber and frequency of the perturbations respectively (ao,6o = const.), we finally obtain the dispersion 
relation 

{n-2kQf = Q\Q^ + 2gul). (3) 

Noting that the plane wave solution is modullationally stable if and only if Im(?7) = 0, we will limit our considerations 
to the focusing case with ^ = — 1, which may potentially lead to MI (apparently, in that case rj becomes complex). 
In this case, Eq. ([2]) shows that if Q < Qcr = UQ^/2 the plane wave solution becomes modulationally unstable. So, 
we notice that there is an implicit dependence on the size of the domain, L. This is, in particular, because the 
wavenumbers accessible in a domain of size L are bounded from below by Qmin = 27r/I/. Therefore, we infer (as was 
also explained in the introduction) that there exists a minimum L, which we denote 

27r V2ii 

^cr=7^ = , (4) 

such that L < L^v =^ Qmin > Qcr, and therefore the uniform solution is not subject to MI and is thus guaranteed to 
be dynamically stable. 



B. Numerical Experiment I: Expanding Domain, No Trap 

It can be anticipated, based on the above discussion, that upon growing such a stable domain beyond the critical size 
Lcr, unstable wavenumbers may emerge. So, starting with L < Lcr, we perturb the plane wave solution of amplitude 
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FIG. 1: Top panel: Space-time evolution plot of the field intensity (square modulus) \uf for an initial condition with uo — 0.5 
and a perturbation of 0.05 cos(27rx). Bottom left panels: Four spatial profiles at times t = 0.2 (stable), t = 0.4 (marginal), 
t = 20 (unstable) and t = 40 (unstable). Bottom right panels: Four respective Fourier profiles; the smallest wavenumber 
present in the solution crosses the critical threshold exactly when the domain size equals the critical value, Lcr = 2v^7r. With 
time, more subcritical wavenumbers appear and are accordingly amplified. 



uq = 0.5 with a stable wavenumber {Q = 27r) and numerically implement the solution of Eq. ([T]) in a growing domain 
{-l{t), l{t)), such that l{t) = 0.5 + lOt (here L = 2/(0) = 1). The resulting evolution is illustrated in Fig. [H The top 
panel shows the space-time evolution of the intensity (the solution's square modulus |i^P), illustrating the emergence 
of a "soliton sprinkler", with an increasing number of solitary waves appearing, as the domain size increases. The 
bottom left four panels show snapshots of the time evolution of the intensity, clearly indicating the growth of the 
perturbation and the eventual formation of multiple solitary wave peaks. Even more importantly, the bottom right set 
of panels elucidates the evolution in Fourier space, with the principal (nonzero) wavenumber initially being above the 
critical wavenumber for the instability (shown by the vertical line). As time evolves however, L increases, hence Qmin 
decreases, and upon the crossing of the critical point we observe growth of the relevant perturbation wavenumbers 
and hence the concomitant generation of the localized peaks in real space. 

We note that in this numerical experiment, as well as in those that will be presented in the following subsections, 
we have used a uniform space mesh of size Ax = 10~^ and a time step of h = 10~^. The methods employed are 
centered, second-order finite difference for spatial discretization and fourth-order Runge-Kutta for time integration. 
We use periodic boundary conditions in order to effectively capture the frequency decomposition. 

In the present implementation, rather than rescaling space as was done in some of the earlier works in the dissipative 
context, we maintain the same mesh-size, while actually increasing the domain with extra nodes at the endpoints of 
the domain. Notice that the above figure (Fig. [T]) doesn't depict the entire final domain, which is approximately 
(—600,600), but merely the relevant activity in the center. 
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C. Numerical Experiment II: Contracting Domain, No Trap 

We now turn to a less orthodox manifestation of MI, which, however, illustrates an important additional dependence 
of the instability's critical point, namely the dependence on the amplitude. If, instead of increasing the domain size, 
we decrease the domain size while increasing background intensity level, so as to maintain the conservation of the 
norm, we also observe a manifestation of ML 

In particular, let us consider a plane wave of amplitude Ui in the finite interval L^^\ Then, decreasing the interval 
to < while conserving the squared norm / \u\^dx^ we obtain a plane wave of amplitude 1^2, which is 

determined by the equation \ui\^L^^^ = |ii2p^*^^^- To this end, utilizing Eq. (|3]), we readily obtain 

LiV=Lil^^<Lil\ (5) 

The above result shows that Qcr^ > Qcr\ which implies that such a shrinking of the domain (while maintaining 
conservation of the norm) has the potential to result in instability of previously stable modes. In fact, for any 

Z, we have a corresponding Qcr = = ui \J 2L(^)/Z. So, as the domain size decreases, the critical wavenumber 
increases, and eventually we expect it to surpass a wavenumber excited in the solution, resulting in MI. This is clearly 
shown in Fig. O which is similar to Fig. (TJ but now for a shrinking domain. The dynamical procedure involves linearly 
decreasing the domain {—l{t)J{t)) with a rate dl/dt = —K and preserving the norm of the modulus numerically 
via the rescaling transformation 

/ \u\^dx= \u\^dx. (6) 

J-lit) J-lo 

We observe that in Fourier space it is now the critical point that is shifting to the right, eventually making the 
configuration unstable, through its crossing of the perturbation wavenumber. The resulting amplification induced by 
the instability is observed both in the real space formation of large amplitude structures and in the Fourier space 
growth of the relevant peaks. A noteworthy feature of the top panel of Fig. [2] is that at t 18 the middle two 
localized structures undergo a nearly elastic collision which underscores their solitary wave character. 



D. Numerical Experiment III: Expanding Domain, Parabolic Trap 

We now turn our attention to the more general inhomogeneous case and consider the full GP equation ([1]) incorpo- 
rating the external potential. The latter is known to be particularly relevant to the confinement of the condensate, and 
assumes typically the harmonic form Ve^t{x) = (l/2)l]^x^, as, e.g., in the case of a magnetic or an optical harmonic 
trap of normalized strength Q pT]. 

Our motivation for considering this example is the following. For Eq. ([1]) with the focusing nonlinearity and 
the parabolic potential, we can obtain the ground state of the system starting, e.g., from the simple linear limit of 
^ = 0. In that case, Eq. ([T]) is actually the usual linear Schrodinger equation, which possesses the exact solutions 
u{x,t) = exp{—i/jjt)ilj{x), where iIj{x) = exp{— Qx'^ / 2^/2) Hn{x), where Hn{x) is the n-th degree Hermite polynomial 
and /i = (2n + l)Cl/^/2 is the normalized chemical potential. One can then straightforwardly establish (based on 
continuation arguments), that there is a nonlinear continuation of this linear state bifurcating from the limit of 
vanishing norm (i.e., vanishing number of atoms). This branch of solutions can be obtained with a fixed point 
algorithm such as a Newton method. There is a unique such continuation for every n; see, e.g., pjj. Here we will 
consider the ground state, emanating from the solution with n = 0. For large amplitude, this solution approaches the 
well-known soliton solution. This is natural since here the nonlinearity tends to focus the solution towards the center 
of the harmonic trap, where Vext{x) and, hence, we revert to the solitonic solution of the untrapped system. 

However, for small amplitudes, this ground state branch represents a more extended, nearly-linear solution (i.e., 
a weakly nonlinear generalization of the linear ground state in the presence of the trap). As such, its width is still 
critically determined by the trap strength Q. Therefore, by changing this parameter we may impose on such a nearly- 
linear solution an effective change in the size of the domain accessible to it. In so doing, we can produce an effectively 
expanding domain, by virtue of decreasing or an effectively contracting domain, by means of increasing Q. These 
are, respectively, the types of numerical experiments that we will present in this and in the following subsections. 

In order to modify the harmonic trap strength, we use a time-dependent parameter 1], according to 



n{t) = 1^1 + {^2 - ^i)(l + tanh[(t - to)/r])/2. 



(7) 
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FIG. 2: Same as Figure [T] but for the decreasing domain case with fixed norm. (Top) The emergence of instability in 
a space-time intensity plot with decreasing domain size beginning with /(O) = 50, K — 10, '0o = 0.5, and a perturbation of 

0. 05 cos(27rx/5). Bottom left: Four typical spatial profiles before (upper) and after (bottom) the manifestation of the instability. 
Bottom right: The four corresponding Fourier profiles are shown, illustrating the emergence of the instability. 

1. e., for t <C to, the trap strength is ^ l^i, for t to, it is ^ 1^2 (and hence we can expand or shrink the domain, 
depending on whether < l^i or 0.2 > ^i, respectively). 

Starting from the case of expanding the domain, we show the relevant results in Figure [3l However, in this version 
of the problem (i.e., with the parabolic trap), when simply expanding the domain the amplitude of the solution will 
decrease (rather than being the same as is the case at the boundaries in Experiment I). However, in order to observe 
the instability, the amplitude of the solution should be preserved. This is achieved by an "injection" of matter in the 
system, so as to preserve the original amplitude. 

This is accomplished in a multiple step process in which we initially integrate until the end of the ramp, at 
approximately t* = to + 3r. Next, we calculate the amplitude loss as a = maxaj(|ii(0)|) — maxa^(|ii(t*)|). Finally we 
rerun the simulation, except this time we augment the right hand side of Eq. ([T]) with a spatially uniform gain rate 
C(t) over the time when the trap opens up. In particular, the considered modified GP equation and the gain term 
assume, respectively, the forms: 

idtu = -d^u^ g\u\'^u^Ve^t{x)u^iC{t)u (8) 
at) = ^sech^[(t-to)/r], (9) 

which may model the employement of a source of BEG atoms, as in the experiment [Ts*]. We then observe that if the 
trap strength is decreased slightly (left set of panels in Fig. [3|), the dynamics is going to be stable, in the sense that 
the distribution of matter will remain unimodal and its width will perform simple small-amplitude oscillations, while 
the Fourier space profile will remain monotonic (quasi-Gaussian, as the Fourier transform of a near-Gaussian profile in 
real space). On the other hand, if the trap strength is decreased significantly (middle and right set of panels in Figure 
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FIG. 3: Top row: Intensity space-time contour plots with Qi = 0.3, = 0.1 (left), = 0.05 (middle), and Q2 = 0.03 (right). 
Middle two rows: The corresponding widths as a function of time (upper panels) and frequencies (lower panels). Bottom two 
rows: Sample spatial profiles (upper panels) and corresponding Fourier profiles (lower panels). 



[3|), then the system becomes "unstable", in the sense that the spatial profile distribution will become multi- modal 
(in fact with more humps developing for smaller frequencies), the oscillations of the matter wave width will become 
more complex (and involve more frequencies) and the Fourier profile will develop nodes, and emergent bands of higher 
wavenumber excitations, whose number is higher the more multi- modal the spatial distribution becomes. 



E. Numerical Experiment IV: Contracting Domain, Parabolic Trap 



Finally, we consider the case where the strength of the harmonic trap is increased for the same type of solutions as 
in Numerical Experiment III. Here, the intensity of the solution is increased as the domain is contracting and, hence 
we can simply allow the dynamics to evolve for different values of Q2 > ^i- What we observe suggests a similar 
scenario of "stability" versus "instability" as in case III. In particular, for small increase of Q (left panels of Fig. 
[4|), the dynamics of the spatial profile remains unimodal (and its Fourier transform near-Gaussian) and the beam 
width oscillates in a near-periodic fashion. On the other hand, for larger trapping strengths (middle and right panels 
of Fig. 2]), the intensity profile involves an increasing number of modes (and its Fourier transform accordingly has 
multiple nodes), and the oscillation of the beam width becomes increasingly irregular, involving a much wider band 
of frequencies. Notice that in this case, similar to case II, no injection of matter is necessary: the instability emerges 
due to the increase in the intensity, necessitated by the conservation of the norm. 
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FIG. 4: Top row: The full solution intensity space-time contour plots for Qi = 0.02, = 0.1 (left), Q2 = 0.3 (middle), and 
^2 = 0.5 (right). Middle two rows: The corresponding widths as a function of time (upper panels) and frequencies (lower 
panels). Bottom two rows: Sample spatial profiles (upper panels) and corresponding Fourier profiles (lower panels). 



III. CONCLUSIONS 



In this work, we have explored how the use of the domain size variations and their effects on instabilities can 
be introduced in a conservative wave setting, by examining the prototypical case of the one-dimensional nonlinear 
Schrodinger equation. We illustrated that this generic model has a modulational instability that can play a role 
similar to the Turing instability of reaction-diffusion systems. In particular, we demonstrated that one can take 
advantage of domain width variability to controllably "excite" the instability and lead to the formation of soliton 
train patterns in the dispersive system. We examined four different manifestations of this feature. Two of them were 
relevant to the uniform, untrapped system and two were associated with the setting of an harmonic trap, relevant to 
Bose-Einstein condensates. In the arguably most standard illustration of the relevant phenomenology (in comparison 
with its dissipative analog), expansion of a domain for a uniform, untrapped solution led to the realization of a "soliton 
sprinkler", emitting more solitons as the domain grew wider. On the other hand, somewhat counter-intuitively, we 
showed that the instability can even emerge for untrapped systems in contracting domains, provided that the "mass" 
(L^ norm) of the profile is preserved. We also demonstrated how the relevant stability or instability notions can 
be re-interpreted in the more realistic, yet not spatially uniform, context of harmonically confined condensates. In 
particular, in the "stable" case, the distribution remained unimodal featuring near-periodic oscillations of its width, 
while in the "unstable" case, multiple humps emerged in the dynamics of the intensity and complex oscillations in its 
moments. This was illustrated both for the expanding domain case of decreasing the parabolic trap strength, as well 
as for the contracting domain case of increasing it. 

It would certainly be relevant to implement similar ideas in higher-dimensional settings, where there also exists the 



8 



additional complication of focusing and wave collapse [8] . It would be especially interesting to elucidate the interplay 
of the different mechanisms, with the modulational instability inducing the emergence of localized waveforms, which, 
in turn, may be subject to collapse. It may also be intriguing to implement such ideas in discrete systems, given 
their recent extensive experimental relevance [1^, Ts] and some of the unique features that they possess such as the 
instability being accessible even for defocusing nonlinearities. Such studies are currently in progress and will be 
reported in future publications. 
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